//
// File: timeDomainHRV.cpp
//
// MATLAB Coder version            : 5.4
// C/C++ source code generated on  : 23-Apr-2025 16:26:18
//

// Include Files
#include "timeDomainHRV.h"
#include "blockedSummation.h"
#include "combineVectorElements.h"
#include "diff.h"
#include "freqDomainHRV_data.h"
#include "freqDomainHRV_initialize.h"
#include "freqDomainHRV_types.h"
#include "hist.h"
#include "linspace.h"
#include "median.h"
#include "minOrMax.h"
#include "rt_nonfinite.h"
#include "std.h"
#include "coder_array.h"
#include <cmath>
#include <cstring>
#include <string.h>

// Function Definitions
//
// timeDomainHRV: calculates time-domain hrv of ibi interval series
//  ibi = 2dim ibi array
//  win = window size to use for sdnni (s)
//  xx = value to use for NNx and pNNx (ms)
//
// Arguments    : const coder::array<double, 2U> &ibi
//                double win
//                double xx
//                struct8_T *output
// Return Type  : void
//
void timeDomainHRV(const coder::array<double, 2U> &ibi, double win, double xx,
                   struct8_T *output)
{
  static const double c[32]{0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
                            1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
                            1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
                            1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.5};
  coder::array<double, 2U> b_peaki_data;
  coder::array<double, 2U> d;
  coder::array<double, 2U> r;
  coder::array<double, 1U> b_ibi;
  coder::array<double, 1U> c_ibi;
  coder::array<double, 1U> tmp;
  coder::array<int, 1U> b_i;
  coder::array<boolean_T, 1U> x;
  double n[32];
  double peaki_data[32];
  double q[32];
  double xout[32];
  double delta1;
  double pi;
  double t_tmp;
  unsigned int a;
  int b_i1;
  int b_i2;
  int b_loop_ub;
  int i;
  int i1;
  int i2;
  int i3;
  int loop_ub;
  int loop_ub_tmp;
  int nx;
  signed char i_data[32];
  boolean_T exitg1;
  if (!isInitialized_freqDomainHRV) {
    freqDomainHRV_initialize();
  }
  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  //  Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
  //
  //  This file is part of HRVAS
  //
  //  HRVAS is free software: you can redistribute it and/or modify
  //  it under the terms of the GNU General Public License as published by
  //  the Free Software Foundation, either version 3 of the License, or
  //  (at your option) any later version.
  //
  //  HRVAS is distributed in the hope that it will be useful,
  //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  //  GNU General Public License for more details.
  //
  //  You should have received a copy of the GNU General Public License
  //  along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  // check inputs
  loop_ub = ibi.size(0);
  b_ibi.set_size(ibi.size(0));
  for (i = 0; i < loop_ub; i++) {
    b_ibi[i] = ibi[i + ibi.size(0)] * 1000.0;
  }
  // convert ibi to ms
  // assumes ibi units are seconds
  //      if abs(range(ibi))<50 %assume ibi units are seconds
  //              ibi=ibi.*1000; %convert ibi to ms
  //      end
  //      if abs(range(diff(t)))>50 %assume time unites are ms
  //          t=t./1000; %convert time to s
  //      end
  //      if t<1000 %assume win units are (s)
  //          t=t*1000; %convert to (ms)
  //      end
  // calculate and round to nearest 1 decimal point
  output->max = std::round(coder::internal::maximum(b_ibi) * 10.0) / 10.0;
  output->min = std::round(coder::internal::minimum(b_ibi) * 10.0) / 10.0;
  output->mean = std::round(coder::blockedSummation(b_ibi, b_ibi.size(0)) /
                            static_cast<double>(b_ibi.size(0)) * 10.0) /
                 10.0;
  output->median = std::round(coder::median(b_ibi) * 10.0) / 10.0;
  output->SDNN = std::round(coder::b_std(b_ibi) * 10.0) / 10.0;
  t_tmp = win * 1000.0;
  // SDNNi: SDNN index is the mean of all the standard deviations of
  // NN (normal RR) intervals for all windows of lenght t.
  a = 0U;
  i1 = 0;
  delta1 = coder::blockedSummation(b_ibi, b_ibi.size(0));
  loop_ub_tmp = static_cast<int>(std::ceil(delta1 / t_tmp));
  tmp.set_size(loop_ub_tmp);
  for (i = 0; i < loop_ub_tmp; i++) {
    tmp[i] = 0.0;
  }
  i = b_ibi.size(0);
  for (i2 = 0; i2 < i; i2++) {
    if (i1 > i2) {
      b_i1 = 0;
      b_i2 = -1;
    } else {
      b_i1 = i1;
      b_i2 = i2;
    }
    loop_ub = b_i2 - b_i1;
    c_ibi.set_size(loop_ub + 1);
    for (b_i2 = 0; b_i2 <= loop_ub; b_i2++) {
      c_ibi[b_i2] = b_ibi[b_i1 + b_i2];
    }
    if (coder::blockedSummation(c_ibi, loop_ub + 1) >= t_tmp) {
      a++;
      if (i1 > i2) {
        i1 = 0;
        b_i1 = -1;
      } else {
        b_i1 = i2;
      }
      loop_ub = b_i1 - i1;
      c_ibi.set_size(loop_ub + 1);
      for (b_i1 = 0; b_i1 <= loop_ub; b_i1++) {
        c_ibi[b_i1] = b_ibi[i1 + b_i1];
      }
      tmp[static_cast<int>(a) - 1] = coder::b_std(c_ibi);
      i1 = i2;
    }
  }
  output->SDANN = std::round(coder::blockedSummation(tmp, tmp.size(0)) /
                             static_cast<double>(tmp.size(0)) * 10.0) /
                  10.0;
  // pNNx: percentage of successive/adjacent NN intervals differing by x (ms) or
  // more
  coder::diff(b_ibi, tmp);
  nx = tmp.size(0);
  c_ibi.set_size(tmp.size(0));
  for (loop_ub = 0; loop_ub < nx; loop_ub++) {
    c_ibi[loop_ub] = std::abs(tmp[loop_ub]);
  }
  // successive ibi diffs (ms)
  x.set_size(c_ibi.size(0));
  loop_ub = c_ibi.size(0);
  for (i = 0; i < loop_ub; i++) {
    x[i] = (c_ibi[i] > xx);
  }
  i1 = x.size(0);
  if (x.size(0) == 0) {
    i2 = 0;
  } else {
    i2 = x[0];
    for (loop_ub = 2; loop_ub <= i1; loop_ub++) {
      i2 += x[loop_ub - 1];
    }
  }
  output->NNx = static_cast<double>(i2) * 10.0 / 10.0;
  i1 = x.size(0);
  if (x.size(0) == 0) {
    i2 = 0;
  } else {
    i2 = x[0];
    for (loop_ub = 2; loop_ub <= i1; loop_ub++) {
      i2 += x[loop_ub - 1];
    }
  }
  output->pNNx = std::round(static_cast<double>(i2) /
                            static_cast<double>(c_ibi.size(0)) * 100.0 * 10.0) /
                 10.0;
  // RMSSD: root mean square of successive RR differences
  // successive ibi diffs
  tmp.set_size(c_ibi.size(0));
  loop_ub = c_ibi.size(0);
  for (i = 0; i < loop_ub; i++) {
    delta1 = c_ibi[i];
    tmp[i] = delta1 * delta1;
  }
  output->RMSSD =
      std::round(std::sqrt(coder::blockedSummation(tmp, tmp.size(0)) /
                           static_cast<double>(c_ibi.size(0))) *
                 10.0) /
      10.0;
  // SDANN: SDANN index is the std of all the mean NN intervals from each
  // segment of lenght t.
  a = 0U;
  i1 = 0;
  tmp.set_size(loop_ub_tmp);
  for (i = 0; i < loop_ub_tmp; i++) {
    tmp[i] = 0.0;
  }
  i = b_ibi.size(0);
  for (i2 = 0; i2 < i; i2++) {
    if (i1 > i2) {
      b_i1 = 0;
      b_i2 = -1;
    } else {
      b_i1 = i1;
      b_i2 = i2;
    }
    loop_ub = b_i2 - b_i1;
    c_ibi.set_size(loop_ub + 1);
    for (b_i2 = 0; b_i2 <= loop_ub; b_i2++) {
      c_ibi[b_i2] = b_ibi[b_i1 + b_i2];
    }
    if (coder::blockedSummation(c_ibi, loop_ub + 1) >= t_tmp) {
      a++;
      if (i1 > i2) {
        i1 = 0;
        b_i1 = -1;
      } else {
        b_i1 = i2;
      }
      loop_ub = b_i1 - i1;
      c_ibi.set_size(loop_ub + 1);
      for (b_i1 = 0; b_i1 <= loop_ub; b_i1++) {
        c_ibi[b_i1] = b_ibi[i1 + b_i1];
      }
      tmp[static_cast<int>(a) - 1] =
          coder::blockedSummation(c_ibi, loop_ub + 1) /
          (static_cast<double>(loop_ub) + 1.0);
      i1 = i2;
    }
  }
  output->SDNNi = std::round(coder::b_std(tmp) * 10.0) / 10.0;
  // heart rate
  tmp.set_size(b_ibi.size(0));
  loop_ub = b_ibi.size(0);
  for (i = 0; i < loop_ub; i++) {
    tmp[i] = 60.0 / (b_ibi[i] / 1000.0);
  }
  output->meanHR = std::round(coder::blockedSummation(tmp, tmp.size(0)) /
                              static_cast<double>(tmp.size(0)) * 10.0) /
                   10.0;
  output->sdHR = std::round(coder::b_std(tmp) * 10.0) / 10.0;
  // GEOMETRIC HRV
  // calculate number of bins to use in histogram
  // 1/128 seconds. Reference: (1996) Heart rate variability: standards of
  // measurement, physiological interpretation and clinical use. temp hrvti: HRV
  // triangular index calculate samples in bin (n) and x location of bins (xout)
  coder::hist(b_ibi, n, xout);
  // hrv ti
  output->HRVTi = std::round(static_cast<double>(b_ibi.size(0)) /
                             coder::internal::maximum(n) * 10.0) /
                  10.0;
  // tinn: triangular interpolation of NN interval histogram
  // Reference: Standards of Measurement, Physiological Interpretation, and
  // Clinical Use
  //            Circulation. 1996; 93(5):1043-1065.
  // calculate histogram of ibi using nbin bins
  coder::hist(b_ibi, n, xout);
  delta1 = coder::internal::maximum(n);
  i1 = 0;
  i2 = 0;
  exitg1 = false;
  while ((!exitg1) && (i2 < 32)) {
    if (n[i2] == delta1) {
      i1++;
      i_data[i1 - 1] = static_cast<signed char>(i2 + 1);
      if (i1 >= 32) {
        exitg1 = true;
      } else {
        i2++;
      }
    } else {
      i2++;
    }
  }
  if (i1 < 1) {
    i1 = 0;
  }
  for (i = 0; i < i1; i++) {
    peaki_data[i] = i_data[i];
  }
  if (i1 > 1) {
    b_peaki_data.set(&peaki_data[0], 1, i1);
    delta1 = std::round(coder::combineVectorElements(b_peaki_data) /
                        static_cast<double>(i1));
    peaki_data[0] = delta1;
  }
  a = 1U;
  pi = peaki_data[0];
  i = static_cast<int>((pi - 1.0) * (32.0 - pi));
  d.set_size(i, 3);
  loop_ub = i * 3;
  for (i = 0; i < loop_ub; i++) {
    d[i] = 0.0;
  }
  i = static_cast<int>(-((-1.0 - (pi - 1.0)) + 1.0));
  for (loop_ub_tmp = 0; loop_ub_tmp < i; loop_ub_tmp++) {
    double m;
    m = (pi - 1.0) + -static_cast<double>(loop_ub_tmp);
    b_i1 = static_cast<int>((1.0 - (pi + 1.0)) + 32.0);
    if (static_cast<int>((1.0 - (pi + 1.0)) + 32.0) - 1 >= 0) {
      if (m < 1.0) {
        b_loop_ub = 0;
      } else {
        b_loop_ub = static_cast<int>(m);
      }
      if (m > pi) {
        i3 = 1;
      } else {
        i3 = static_cast<int>(m);
      }
    }
    for (int b_n{0}; b_n < b_i1; b_n++) {
      double y[32];
      double c_n;
      c_n = (pi + 1.0) + static_cast<double>(b_n);
      // define triangle that fits the histogram
      std::memset(&q[0], 0, 32U * sizeof(double));
      if (b_loop_ub - 1 >= 0) {
        std::memset(&q[0], 0, b_loop_ub * sizeof(double));
      }
      if (c_n > 32.0) {
        b_i2 = -1;
        nx = -1;
      } else {
        b_i2 = static_cast<int>(c_n) - 2;
        nx = 31;
      }
      loop_ub = nx - b_i2;
      if (loop_ub - 1 >= 0) {
        std::memset(&q[b_i2 + 1], 0,
                    ((loop_ub + b_i2) - b_i2) * sizeof(double));
      }
      t_tmp = n[static_cast<int>(pi) - 1];
      coder::linspace(t_tmp, (pi - m) + 1.0, r);
      loop_ub = r.size(1);
      for (b_i2 = 0; b_i2 < loop_ub; b_i2++) {
        q[(i3 + b_i2) - 1] = r[b_i2];
      }
      if (pi > c_n) {
        b_i2 = 1;
      } else {
        b_i2 = static_cast<int>(pi);
      }
      delta1 = (c_n - pi) + 1.0;
      if (!(delta1 >= 0.0)) {
        r.set_size(1, 0);
      } else {
        r.set_size(1, static_cast<int>(delta1));
        if (static_cast<int>(delta1) >= 1) {
          nx = static_cast<int>(delta1) - 1;
          r[static_cast<int>(delta1) - 1] = 0.0;
          if (r.size(1) >= 2) {
            r[0] = t_tmp;
            if (r.size(1) >= 3) {
              if ((t_tmp == -0.0) && (static_cast<int>(delta1) > 2)) {
                for (loop_ub = 2; loop_ub <= nx; loop_ub++) {
                  r[loop_ub - 1] =
                      (static_cast<double>((loop_ub << 1) -
                                           static_cast<int>(delta1)) -
                       1.0) *
                      0.0;
                }
                if ((static_cast<int>(delta1) & 1) == 1) {
                  r[static_cast<int>(delta1) >> 1] = 0.0;
                }
              } else {
                delta1 = (0.0 - t_tmp) / (static_cast<double>(r.size(1)) - 1.0);
                nx = r.size(1);
                for (loop_ub = 0; loop_ub <= nx - 3; loop_ub++) {
                  r[loop_ub + 1] =
                      t_tmp + (static_cast<double>(loop_ub) + 1.0) * delta1;
                }
              }
            }
          }
        }
      }
      loop_ub = r.size(1);
      for (nx = 0; nx < loop_ub; nx++) {
        q[(b_i2 + nx) - 1] = r[nx];
      }
      // integrate squared difference
      for (loop_ub = 0; loop_ub < 32; loop_ub++) {
        t_tmp = n[loop_ub] - q[loop_ub];
        q[loop_ub] = t_tmp;
        y[loop_ub] = t_tmp * t_tmp;
      }
      t_tmp = 0.0;
      for (nx = 0; nx < 32; nx++) {
        i2 = nx + 1;
        delta1 = c[nx];
        for (i1 = i2; i1 <= i2; i1++) {
          t_tmp += y[i1 - 1] * delta1;
        }
      }
      i1 = static_cast<int>(a + b_n) - 1;
      d[i1] = t_tmp;
      d[i1 + d.size(0)] = m;
      d[i1 + d.size(0) * 2] = c_n;
      // plot(D); hold on; plot(q,'r'); hold off;
      // title(['d^2 = ' num2str(d(i,1))])
    }
    a += static_cast<int>((1.0 - (pi + 1.0)) + 32.0);
  }
  // find where minimum square diff occured
  loop_ub = d.size(0);
  c_ibi.set_size(d.size(0));
  for (i = 0; i < loop_ub; i++) {
    c_ibi[i] = d[i];
  }
  delta1 = coder::internal::minimum(c_ibi);
  loop_ub = d.size(0);
  x.set_size(d.size(0));
  for (i = 0; i < loop_ub; i++) {
    x[i] = (d[i] == delta1);
  }
  nx = x.size(0);
  i1 = 0;
  b_i.set_size(x.size(0));
  i2 = 0;
  exitg1 = false;
  while ((!exitg1) && (i2 <= nx - 1)) {
    if (x[i2]) {
      i1++;
      b_i[i1 - 1] = i2 + 1;
      if (i1 >= nx) {
        exitg1 = true;
      } else {
        i2++;
      }
    } else {
      i2++;
    }
  }
  if (x.size(0) == 1) {
    if (i1 == 0) {
      b_i.set_size(0);
    }
  } else {
    if (i1 < 1) {
      i1 = 0;
    }
    b_i.set_size(i1);
  }
  tmp.set_size(b_i.size(0));
  loop_ub = b_i.size(0);
  for (i = 0; i < loop_ub; i++) {
    tmp[i] = b_i[i];
  }
  // make sure there is only one choise
  // calculate TINN in (ms)
  // plot
  //      X=xout(peaki); M=xout(m); N=xout(n); Y=nout(peaki);
  //      figure;
  //      hist(ibi,nbin)
  //      xlimits=get(gca,'xlim');
  //      hold on;
  //      plot(xout,nout,'k')
  //      line([M X N M],[0 Y 0 0],'color','r','linewidth',1.5,'LineStyle','--')
  //      line([X X],[0 1000],'LineStyle','-.','color','k')
  //      line([0 2000],[Y Y],'LineStyle','-.','color','k')
  //      colormap white
  //
  //      xlabel('IBI (ms)');
  //      ylabel('Number of IBI')
  //      legend({'Histogram','D(t)','q(t)'})
  //      set(gca,'xtick',[xout(m) xout(peaki)
  //      xout(n)],'xticklabel',{'N','X','M'}, ...
  //          'ytick',Y,'yticklabel','Y')
  //      set(gca,'xlim',xlimits);
  output->TINN =
      std::round(
          std::abs(xout[static_cast<int>(
                            d[(static_cast<int>(tmp[0]) + d.size(0) * 2) - 1]) -
                        1] -
                   xout[static_cast<int>(
                            d[(static_cast<int>(tmp[0]) + d.size(0)) - 1]) -
                        1]) *
          10.0) /
      10.0;
}

//
// File trailer for timeDomainHRV.cpp
//
// [EOF]
//
